Simon Frost (@sdwfrost), 2020-05-20
One high-level representation of the SIR model is as a reaction network, borrowed from systems biology. ModelingToolkit allows us to convert this representation to ODEs, SDEs, and jump processes. This example is a slightly tweaked version of one in the ModelingToolkit documentation, using the population size as a derived variable in the rates for the transitions.
using DifferentialEquations using ModelingToolkit using OrdinaryDiffEq using StochasticDiffEq using DiffEqJump using Random using Plots
@parameters t β c γ @variables S(t) I(t) R(t) N=S+I+R # This is recognized as a derived variable rxs = [Reaction((β*c)/N, [S,I], [I], [1,1], [2]) Reaction(γ, [I], [R])]
2-element Array{ModelingToolkit.Reaction{Any,Int64},1}:
ModelingToolkit.Reaction{Any,Int64}(c*β*(((I(t)) + (R(t)) + (S(t)))^-1), S
ymbolicUtils.Term{Real}[S(t), I(t)], SymbolicUtils.Term{Real}[I(t)], [1, 1]
, [2], Pair{Any,Int64}[S(t) => -1, I(t) => 1], false)
ModelingToolkit.Reaction{Any,Int64}(γ, SymbolicUtils.Term{Real}[I(t)], Sym
bolicUtils.Term{Real}[R(t)], [1], [1], Pair{Any,Int64}[I(t) => -1, R(t) =>
1], false)
rs = ReactionSystem(rxs, t, [S,I,R], [β,c,γ])
Model ##ReactionSystem#256 with 2 equations States (3): S(t) I(t) R(t) Parameters (3): β c γ
We set the timespan for simulations, tspan, initial conditions, u0, and parameter values, p (which are unpacked above as [β,γ]).
tmax = 40.0 tspan = (0.0,tmax);
In ModelingToolkit, the initial values are defined by a dictionary.
u0 = [S => 990.0, I => 10.0, R => 0.0];
Similarly, the parameter values are defined by a dictionary.
p = [β=>0.05, c=>10.0, γ=>0.25];
Random.seed!(1234);
odesys = convert(ODESystem, rs) oprob = ODEProblem(odesys, u0, tspan, p) osol = solve(oprob, Tsit5()) plot(osol)
sdesys = convert(SDESystem, rs) sprob = SDEProblem(sdesys, u0, tspan, p) ssol = solve(sprob, LambaEM()) plot(ssol)
To convert to a jump process, we need to set the initial conditions to Int rather than Float.
jumpsys = convert(JumpSystem, rs) u0i = [S => 990, I => 10, R => 0] dprob = DiscreteProblem(jumpsys, u0i, tspan, p) jprob = JumpProblem(jumpsys, dprob, Direct()) jsol = solve(jprob, SSAStepper()) plot(jsol)